source('../env.R')
community_data = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics_using_relative_abundance.csv'))
head(community_data)
colnames(community_data)
[1] "mntd_standard" "mntd_normal" "mntd_actual" "mass_fdiv_standard"
[5] "mass_fdiv_normal" "mass_fdiv_actual" "beak_width_fdiv_standard" "beak_width_fdiv_normal"
[9] "beak_width_fdiv_actual" "hwi_fdiv_standard" "hwi_fdiv_normal" "hwi_fdiv_actual"
[13] "city_id" "urban_pool_size"
min(community_data$mntd_standard)
[1] -2.338219
max(community_data$mntd_standard)
[1] 2.34743
min(community_data$beak_width_fdiv_standard)
[1] -2.685152
max(community_data$beak_width_fdiv_standard)
[1] 1.931681
min(community_data$hwi_fdiv_standard)
[1] -2.254896
max(community_data$hwi_fdiv_standard)
[1] 2.315362
min(community_data$mass_fdiv_standard)
[1] -2.351968
max(community_data$mass_fdiv_standard)
[1] 2.114779
min(community_data$mntd_normal)
[1] 0
max(community_data$mntd_normal)
[1] 1
min(community_data$beak_width_fdiv_normal)
[1] 0
max(community_data$beak_width_fdiv_normal)
[1] 1
min(community_data$hwi_fdiv_normal)
[1] 0
max(community_data$hwi_fdiv_normal)
[1] 1
min(community_data$mass_fdiv_normal)
[1] 0
max(community_data$mass_fdiv_normal)
[1] 1
Join on realms
city_to_realm = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))
community_data_with_realm = left_join(community_data, city_to_realm)
Cities as points
city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp'))) %>% left_join(community_data_with_realm)
city_points_coords = st_coordinates(city_points)
city_points$latitude = city_points_coords[,1]
city_points$longitude = city_points_coords[,2]
world_map = read_country_boundaries()
Load community data, and create long format version
communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))
communities
community_summary = communities %>% group_by(city_id) %>% summarise(regional_pool_size = n(), urban_pool_size = sum(relative_abundance_proxy > 0))
community_summary
Load trait data
traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_ebird.csv'))
head(traits)
Load spatial var
spatial_var = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'spatial_var.csv')) %>% filter(city_id %in% community_summary$city_id)
spatial_var
test_required_values = function(name, df) {
cat(paste(
test_value_wilcox(paste(name, 'Std: MNTD'), df$mntd_standard),
test_value_wilcox(paste(name, 'Std: Beak Width FDiv'), df$beak_width_fdiv_standard),
test_value_wilcox(paste(name, 'Std: HWI FDiv'), df$hwi_fdiv_standard),
test_value_wilcox(paste(name, 'Std: Mass FDiv'), df$mass_fdiv_standard),
test_value_wilcox(paste(name, 'Norm: MNTD'), df$mntd_normal, mu = 0.5),
test_value_wilcox(paste(name, 'Norm: Beak Width FDiv'), df$beak_width_fdiv_normal, mu = 0.5),
test_value_wilcox(paste(name, 'Norm: HWI FDiv'), df$hwi_fdiv_normal, mu = 0.5),
test_value_wilcox(paste(name, 'Norm: Mass FDiv'), df$mass_fdiv_normal, mu = 0.5),
paste('N', nrow(df)),
sep = "\n"))
}
test_required_values('Global', community_data_with_realm)
Global Std: MNTD median -0.36 ***
Global Std: Beak Width FDiv median 0.03
Global Std: HWI FDiv median 0.39 **
Global Std: Mass FDiv median 0.29 ***
Global Norm: MNTD median 0.38 ***
Global Norm: Beak Width FDiv median 0.59 ***
Global Norm: HWI FDiv median 0.64 ***
Global Norm: Mass FDiv median 0.64 ***
N 308
unique(community_data_with_realm$core_realm)
[1] "Nearctic" "Neotropic" "Palearctic" "Afrotropic" "Indomalayan" "Australasia"
test_required_values('Nearctic', community_data_with_realm[community_data_with_realm$core_realm == 'Nearctic',])
Nearctic Std: MNTD median 0.67 *
Nearctic Std: Beak Width FDiv median 0.29
Nearctic Std: HWI FDiv median -0.8 ***
Nearctic Std: Mass FDiv median -0.26
Nearctic Norm: MNTD median 0.75
Nearctic Norm: Beak Width FDiv median 0.62
Nearctic Norm: HWI FDiv median 0.23 **
Nearctic Norm: Mass FDiv median 0.48
N 46
test_required_values('Neotropic', community_data_with_realm[community_data_with_realm$core_realm == 'Neotropic',])
Neotropic Std: MNTD median 0.03
Neotropic Std: Beak Width FDiv median -0.44 ***
Neotropic Std: HWI FDiv median -0.31
Neotropic Std: Mass FDiv median 0.33 *
Neotropic Norm: MNTD median 0.45
Neotropic Norm: Beak Width FDiv median 0.46
Neotropic Norm: HWI FDiv median 0.5
Neotropic Norm: Mass FDiv median 0.65 ***
N 64
test_required_values('Palearctic', community_data_with_realm[community_data_with_realm$core_realm == 'Palearctic',])
Palearctic Std: MNTD median 0.13
Palearctic Std: Beak Width FDiv median 1.25 ***
Palearctic Std: HWI FDiv median -0.39
Palearctic Std: Mass FDiv median 0.01
Palearctic Norm: MNTD median 0.57 *
Palearctic Norm: Beak Width FDiv median 0.93 ***
Palearctic Norm: HWI FDiv median 0.38
Palearctic Norm: Mass FDiv median 0.55
N 72
test_required_values('Afrotropic', community_data_with_realm[community_data_with_realm$core_realm == 'Afrotropic',])
Afrotropic Std: MNTD median -1.27 *
Afrotropic Std: Beak Width FDiv median -0.5
Afrotropic Std: HWI FDiv median 0.17
Afrotropic Std: Mass FDiv median -0.95
Afrotropic Norm: MNTD median 0.07 *
Afrotropic Norm: Beak Width FDiv median 0.42
Afrotropic Norm: HWI FDiv median 0.54
Afrotropic Norm: Mass FDiv median 0.34
N 9
test_required_values('Indomalayan', community_data_with_realm[community_data_with_realm$core_realm == 'Indomalayan',])
Indomalayan Std: MNTD median -0.63 ***
Indomalayan Std: Beak Width FDiv median -0.7 ***
Indomalayan Std: HWI FDiv median 1.11 ***
Indomalayan Std: Mass FDiv median 0.84 ***
Indomalayan Norm: MNTD median 0.21 ***
Indomalayan Norm: Beak Width FDiv median 0.46
Indomalayan Norm: HWI FDiv median 0.88 ***
Indomalayan Norm: Mass FDiv median 0.8 ***
N 111
test_required_values('Australasia', community_data_with_realm[community_data_with_realm$core_realm == 'Australasia',])
Australasia Std: MNTD median -1.39
Australasia Std: Beak Width FDiv median -0.74
Australasia Std: HWI FDiv median 0.78
Australasia Std: Mass FDiv median -0.94
Australasia Norm: MNTD median 0.09
Australasia Norm: Beak Width FDiv median 0.46
Australasia Norm: HWI FDiv median 0.79
Australasia Norm: Mass FDiv median 0.55
N 6
print('Standard')
[1] "Standard"
kruskal.test(mntd_standard~core_realm, data = community_data_with_realm)
Kruskal-Wallis rank sum test
data: mntd_standard by core_realm
Kruskal-Wallis chi-squared = 102.11, df = 5, p-value < 0.00000000000000022
pairwise.wilcox.test(community_data_with_realm$mntd_standard, community_data_with_realm$core_realm)
Pairwise comparisons using Wilcoxon rank sum exact test
data: community_data_with_realm$mntd_standard and community_data_with_realm$core_realm
Afrotropic Australasia Indomalayan Nearctic Neotropic
Australasia 1.00000 - - - -
Indomalayan 0.00489 0.01024 - - -
Nearctic 0.0000053220919 0.00018 0.0000000141740 - -
Neotropic 0.00018 0.00170 0.0000000083499 0.21142 -
Palearctic 0.00018 0.00170 0.0000000000087 0.04492 1.00000
P value adjustment method: holm
print('Normal')
[1] "Normal"
kruskal.test(mntd_normal~core_realm, data = community_data_with_realm)
Kruskal-Wallis rank sum test
data: mntd_normal by core_realm
Kruskal-Wallis chi-squared = 122.2, df = 5, p-value < 0.00000000000000022
pairwise.wilcox.test(community_data_with_realm$mntd_normal, community_data_with_realm$core_realm)
Pairwise comparisons using Wilcoxon rank sum exact test
data: community_data_with_realm$mntd_normal and community_data_with_realm$core_realm
Afrotropic Australasia Indomalayan Nearctic Neotropic
Australasia 0.27213 - - - -
Indomalayan 0.00013 0.00267 - - -
Nearctic 0.00276 0.01815 0.0000009331952 - -
Neotropic 0.0000199762545 0.00053 0.0000000000089 0.01638 -
Palearctic 0.0000540321246 0.00147 < 0.0000000000000002 0.10733 0.00947
P value adjustment method: holm
print('Standard')
[1] "Standard"
kruskal.test(beak_width_fdiv_standard~core_realm, data = community_data_with_realm)
Kruskal-Wallis rank sum test
data: beak_width_fdiv_standard by core_realm
Kruskal-Wallis chi-squared = 107.95, df = 5, p-value < 0.00000000000000022
pairwise.wilcox.test(community_data_with_realm$beak_width_fdiv_standard, community_data_with_realm$core_realm)
Pairwise comparisons using Wilcoxon rank sum exact test
data: community_data_with_realm$beak_width_fdiv_standard and community_data_with_realm$core_realm
Afrotropic Australasia Indomalayan Nearctic Neotropic
Australasia 1.00000 - - - -
Indomalayan 1.00000 1.00000 - - -
Nearctic 0.11361 1.00000 0.00123 - -
Neotropic 1.00000 1.00000 1.00000 0.00388 -
Palearctic 0.00026 0.13449 < 0.0000000000000002 0.000001461233727 0.000000000000018
P value adjustment method: holm
print('Normal')
[1] "Normal"
kruskal.test(beak_width_fdiv_standard~core_realm, data = community_data_with_realm)
Kruskal-Wallis rank sum test
data: beak_width_fdiv_standard by core_realm
Kruskal-Wallis chi-squared = 107.95, df = 5, p-value < 0.00000000000000022
pairwise.wilcox.test(community_data_with_realm$beak_width_fdiv_standard, community_data_with_realm$core_realm)
Pairwise comparisons using Wilcoxon rank sum exact test
data: community_data_with_realm$beak_width_fdiv_standard and community_data_with_realm$core_realm
Afrotropic Australasia Indomalayan Nearctic Neotropic
Australasia 1.00000 - - - -
Indomalayan 1.00000 1.00000 - - -
Nearctic 0.11361 1.00000 0.00123 - -
Neotropic 1.00000 1.00000 1.00000 0.00388 -
Palearctic 0.00026 0.13449 < 0.0000000000000002 0.000001461233727 0.000000000000018
P value adjustment method: holm
print('Standard')
[1] "Standard"
kruskal.test(hwi_fdiv_standard~core_realm, data = community_data_with_realm)
Kruskal-Wallis rank sum test
data: hwi_fdiv_standard by core_realm
Kruskal-Wallis chi-squared = 115.62, df = 5, p-value < 0.00000000000000022
pairwise.wilcox.test(community_data_with_realm$hwi_fdiv_standard, community_data_with_realm$core_realm)
Pairwise comparisons using Wilcoxon rank sum exact test
data: community_data_with_realm$hwi_fdiv_standard and community_data_with_realm$core_realm
Afrotropic Australasia Indomalayan Nearctic Neotropic
Australasia 1.0000 - - - -
Indomalayan 0.1471 1.0000 - - -
Nearctic 1.0000 0.0038 0.00000000000000043 - -
Neotropic 1.0000 0.0680 0.00000000000000197 0.0925 -
Palearctic 1.0000 0.3068 0.00000000003605628 0.2673 1.0000
P value adjustment method: holm
print('Normal')
[1] "Normal"
kruskal.test(hwi_fdiv_normal~core_realm, data = community_data_with_realm)
Kruskal-Wallis rank sum test
data: hwi_fdiv_normal by core_realm
Kruskal-Wallis chi-squared = 104.9, df = 5, p-value < 0.00000000000000022
pairwise.wilcox.test(community_data_with_realm$hwi_fdiv_normal, community_data_with_realm$core_realm)
Pairwise comparisons using Wilcoxon rank sum exact test
data: community_data_with_realm$hwi_fdiv_normal and community_data_with_realm$core_realm
Afrotropic Australasia Indomalayan Nearctic Neotropic
Australasia 0.7233 - - - -
Indomalayan 0.1144 1.0000 - - -
Nearctic 0.0881 0.0193 0.0000000000006844 - -
Neotropic 1.0000 0.0204 0.0000000000000014 0.0000964104748761 -
Palearctic 0.9971 0.2748 0.0000000522948450 0.0073 0.9953
P value adjustment method: holm
print('Standard')
[1] "Standard"
kruskal.test(mass_fdiv_standard~core_realm, data = community_data_with_realm)
Kruskal-Wallis rank sum test
data: mass_fdiv_standard by core_realm
Kruskal-Wallis chi-squared = 48.784, df = 5, p-value = 0.000000002457
pairwise.wilcox.test(community_data_with_realm$mass_fdiv_standard, community_data_with_realm$core_realm)
Pairwise comparisons using Wilcoxon rank sum exact test
data: community_data_with_realm$mass_fdiv_standard and community_data_with_realm$core_realm
Afrotropic Australasia Indomalayan Nearctic Neotropic
Australasia 1.0000 - - - -
Indomalayan 0.0033 0.1401 - - -
Nearctic 0.1401 1.0000 0.0023 - -
Neotropic 0.0158 0.2989 1.0000 0.1401 -
Palearctic 0.0873 1.0000 0.00000002 1.0000 0.0658
P value adjustment method: holm
print('Normal')
[1] "Normal"
kruskal.test(mass_fdiv_normal~core_realm, data = community_data_with_realm)
Kruskal-Wallis rank sum test
data: mass_fdiv_normal by core_realm
Kruskal-Wallis chi-squared = 58.138, df = 5, p-value = 0.00000000002946
pairwise.wilcox.test(community_data_with_realm$mass_fdiv_normal, community_data_with_realm$core_realm)
Pairwise comparisons using Wilcoxon rank sum exact test
data: community_data_with_realm$mass_fdiv_normal and community_data_with_realm$core_realm
Afrotropic Australasia Indomalayan Nearctic Neotropic
Australasia 1.000 - - - -
Indomalayan 0.014 0.075 - - -
Nearctic 1.000 1.000 0.016 - -
Neotropic 0.126 1.000 0.014 0.958 -
Palearctic 0.824 1.000 0.0000000000008 1.000 0.014
P value adjustment method: holm
communities %>%
left_join(city_to_realm) %>%
mutate(family = gsub( " .*$", "", ebird_species_name)) %>%
dplyr::select(family, core_realm) %>%
distinct() %>%
arrange(core_realm)
communities %>%
mutate(family = gsub( " .*$", "", ebird_species_name)) %>%
dplyr::select(family) %>%
distinct() %>%
arrange()
of which urban
communities %>%
filter(relative_abundance_proxy > 0) %>%
mutate(family = gsub( " .*$", "", ebird_species_name)) %>%
dplyr::select(family) %>%
distinct() %>%
arrange()
communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))
city_introduced_species = communities %>% group_by(city_id) %>% summarise(number_of_species = n()) %>% left_join(
communities %>% group_by(city_id) %>% filter(origin == 'Introduced') %>% summarise(number_of_introduced_species = n())
) %>% replace_na(list(number_of_introduced_species = 0))
community_data_with_introductions = left_join(community_data, city_introduced_species)
community_data_with_introductions$has_introduced_species = community_data_with_introductions$number_of_introduced_species > 0
community_data_with_introductions
communities %>%
filter(origin == 'Introduced') %>%
dplyr::select(ebird_species_name) %>%
group_by(ebird_species_name) %>%
summarise(total_cities = n()) %>%
arrange(desc(total_cities))
community_data_with_introductions[,c('mntd_standard', 'has_introduced_species')]
community_data_with_introductions %>% group_by(has_introduced_species) %>% summarise(
total_cities = n(),
mean_mntd_std = mean(mntd_standard, na.rm = T),
median_mntd_std = median(mntd_standard, na.rm = T),
sd_mntd_std = sd(mntd_standard, na.rm = T),
mean_mass_fdiv_std = mean(mass_fdiv_standard, na.rm = T),
median_mass_fdiv_std = median(mass_fdiv_standard, na.rm = T),
sd_mass_fdiv_std = sd(mass_fdiv_standard, na.rm = T),
mean_gape_width_fdiv_std = mean(beak_width_fdiv_standard, na.rm = T),
median_gape_width_fdiv_std = median(beak_width_fdiv_standard, na.rm = T),
sd_gape_width_fdiv_std = sd(beak_width_fdiv_standard, na.rm = T),
mean_handwing_index_fdiv_std = mean(hwi_fdiv_standard, na.rm = T),
median_handwing_index_fdiv_std = median(hwi_fdiv_standard, na.rm = T),
sd_handwing_index_fdiv_std = sd(hwi_fdiv_standard, na.rm = T)
)
ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = mntd_standard)) + geom_boxplot()
ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = mntd_normal)) + geom_boxplot()
wilcox.test(mntd_standard ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')
Wilcoxon rank sum test with continuity correction
data: mntd_standard by has_introduced_species
W = 7922, p-value = 0.00001262
alternative hypothesis: true location shift is not equal to 0
There is a significant difference between the response of cities with introduced species (0.53±0.27) and those without (0.47±0.19) (p-value = 0.02).
wilcox.test(mntd_normal ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')
Wilcoxon rank sum test with continuity correction
data: mntd_normal by has_introduced_species
W = 7376.5, p-value = 0.0000003705
alternative hypothesis: true location shift is not equal to 0
ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = mass_fdiv_standard)) + geom_boxplot()
ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = mass_fdiv_normal)) + geom_boxplot()
wilcox.test(mass_fdiv_standard ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')
Wilcoxon rank sum test with continuity correction
data: mass_fdiv_standard by has_introduced_species
W = 15009, p-value = 0.0000007626
alternative hypothesis: true location shift is not equal to 0
There is a significant difference between the response of cities with introduced species (0.57±0.27) and those without (0.73±0.24) (p < 0.0001)
wilcox.test(mass_fdiv_normal ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')
Wilcoxon rank sum test with continuity correction
data: mass_fdiv_normal by has_introduced_species
W = 15204, p-value = 0.0000001989
alternative hypothesis: true location shift is not equal to 0
ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = beak_width_fdiv_standard)) + geom_boxplot()
ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = beak_width_fdiv_normal)) + geom_boxplot()
wilcox.test(beak_width_fdiv_standard ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')
Wilcoxon rank sum test with continuity correction
data: beak_width_fdiv_standard by has_introduced_species
W = 8659, p-value = 0.0006786
alternative hypothesis: true location shift is not equal to 0
There is NOT a significant difference between the response of cities with introduced species (0.61±0.30) and those without (0.56±0.27)
wilcox.test(beak_width_fdiv_normal ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')
Wilcoxon rank sum test with continuity correction
data: beak_width_fdiv_normal by has_introduced_species
W = 9829.5, p-value = 0.06288
alternative hypothesis: true location shift is not equal to 0
ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = hwi_fdiv_standard)) + geom_boxplot()
ggplot(community_data_with_introductions, aes(x = has_introduced_species, y = hwi_fdiv_normal)) + geom_boxplot()
wilcox.test(hwi_fdiv_standard ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')
Wilcoxon rank sum test with continuity correction
data: hwi_fdiv_standard by has_introduced_species
W = 17622, p-value < 0.00000000000000022
alternative hypothesis: true location shift is not equal to 0
There is a significant difference between the response of cities with introduced species (0.49±0.30) and those without (0.79±0.21) (p < 0.0001)
wilcox.test(hwi_fdiv_normal ~ has_introduced_species, community_data_with_introductions, na.action = 'na.omit')
Wilcoxon rank sum test with continuity correction
data: hwi_fdiv_normal by has_introduced_species
W = 17673, p-value < 0.00000000000000022
alternative hypothesis: true location shift is not equal to 0
community_data_with_introductions %>% left_join(city_to_realm) %>%
group_by(core_realm) %>%
summarise(
total_cities = n(),
total_cities_with_introduced = sum(has_introduced_species),
total_cities_without_introduced = n() - sum(has_introduced_species)) %>%
arrange(core_realm)
communities %>%
filter(origin == 'Introduced') %>%
filter(relative_abundance_proxy < 0.1)
communities %>%
group_by(origin) %>%
summarise(average_relative_abundance = mean(relative_abundance_proxy))
communities %>%
group_by(origin) %>%
filter(relative_abundance_proxy > 0) %>%
summarise(average_relative_abundance = mean(relative_abundance_proxy))
communities %>%
group_by(origin) %>%
summarise(average_relative_abundance = mean(relative_abundance_proxy))
geography = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'geography.csv'))
names(geography)
[1] "city_id" "city_avg_ndvi" "city_avg_elevation" "city_avg_temp"
[5] "city_avg_min_monthly_temp" "city_avg_max_monthly_temp" "city_avg_monthly_temp" "city_avg_rainfall"
[9] "city_avg_max_monthly_rainfall" "city_avg_min_monthly_rainfall" "city_avg_soil_moisture" "city_max_elev"
[13] "city_min_elev" "city_elev_range" "region_20km_avg_ndvi" "region_20km_avg_elevation"
[17] "region_20km_avg_soil_moisture" "region_20km_max_elev" "region_20km_min_elev" "region_20km_elev_range"
[21] "region_50km_avg_ndvi" "region_50km_avg_elevation" "region_50km_avg_soil_moisture" "region_50km_max_elev"
[25] "region_50km_min_elev" "region_50km_elev_range"
analysis_data = community_data_with_realm[,c('city_id',
'mntd_standard', 'mass_fdiv_standard', 'beak_width_fdiv_standard', 'hwi_fdiv_standard',
'mntd_normal', 'mass_fdiv_normal', 'beak_width_fdiv_normal', 'hwi_fdiv_normal',
'core_realm')] %>%
left_join(city_points[,c('city_id', 'latitude', 'longitude')]) %>%
left_join(community_data_with_introductions[,c('city_id', 'has_introduced_species')]) %>%
left_join(geography) %>%
left_join(spatial_var)
analysis_data$abs_latitude = abs(analysis_data$latitude)
analysis_data$core_realm = factor(analysis_data$core_realm, levels = c('Palearctic', 'Nearctic', 'Neotropic', 'Afrotropic', 'Indomalayan', 'Australasia', 'Oceania'))
analysis_data$has_introduced_species = factor(analysis_data$has_introduced_species, level = c('FALSE', 'TRUE'), labels = c('No introduced species', 'Introduced species'))
model_data = function(df, dependant_var) {
df[,c(dependant_var, 'core_realm', 'abs_latitude', 'longitude', 'has_introduced_species', 'city_avg_ndvi', 'city_avg_elevation', 'city_avg_temp', 'city_avg_min_monthly_temp', 'city_avg_max_monthly_temp', 'city_avg_monthly_temp', 'city_avg_rainfall', 'city_avg_max_monthly_rainfall', 'city_avg_min_monthly_rainfall', 'city_avg_soil_moisture', 'city_max_elev', 'city_min_elev', 'city_elev_range', 'region_20km_avg_ndvi', 'region_20km_avg_elevation', 'region_20km_avg_soil_moisture', 'region_20km_max_elev', 'region_20km_min_elev', 'region_20km_elev_range', 'region_50km_avg_ndvi', 'region_50km_avg_elevation', 'region_50km_avg_soil_moisture', 'region_50km_max_elev', 'region_50km_min_elev', 'region_50km_elev_range')]
}
model_data(analysis_data, 'mntd_standard')
model_data(analysis_data, 'mntd_normal')
names(analysis_data)
[1] "city_id" "mntd_standard" "mass_fdiv_standard" "beak_width_fdiv_standard"
[5] "hwi_fdiv_standard" "mntd_normal" "mass_fdiv_normal" "beak_width_fdiv_normal"
[9] "hwi_fdiv_normal" "core_realm" "latitude" "longitude"
[13] "geometry" "has_introduced_species" "city_avg_ndvi" "city_avg_elevation"
[17] "city_avg_temp" "city_avg_min_monthly_temp" "city_avg_max_monthly_temp" "city_avg_monthly_temp"
[21] "city_avg_rainfall" "city_avg_max_monthly_rainfall" "city_avg_min_monthly_rainfall" "city_avg_soil_moisture"
[25] "city_max_elev" "city_min_elev" "city_elev_range" "region_20km_avg_ndvi"
[29] "region_20km_avg_elevation" "region_20km_avg_soil_moisture" "region_20km_max_elev" "region_20km_min_elev"
[33] "region_20km_elev_range" "region_50km_avg_ndvi" "region_50km_avg_elevation" "region_50km_avg_soil_moisture"
[37] "region_50km_max_elev" "region_50km_min_elev" "region_50km_elev_range" "NMDS1"
[41] "NMDS2" "abs_latitude"
analysis_data_nmds_coords = analysis_data[,c('NMDS1', 'NMDS2')]
coordinates(analysis_data_nmds_coords) = ~ NMDS1 + NMDS2
analysis_data_nmds_nearneigh <- knearneigh(analysis_data_nmds_coords)
analysis_data_nmds_neighbours <- knn2nb(analysis_data_nmds_nearneigh)
cities_to_realms_nmds = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv')) %>% left_join(analysis_data) %>% filter(!is.na(NMDS1))
unique(cities_to_realms_nmds$core_realm)
[1] "Nearctic" "Neotropic" "Palearctic" "Afrotropic" "Indomalayan" "Australasia"
realm_nmds_neartic_polygon = cities_to_realms_nmds %>% filter(core_realm == 'Nearctic') %>% slice(chull(NMDS1, NMDS2))
realm_nmds_neotropic_polygon = cities_to_realms_nmds %>% filter(core_realm == 'Neotropic') %>% slice(chull(NMDS1, NMDS2))
realm_nmds_palearctic_polygon = cities_to_realms_nmds %>% filter(core_realm == 'Palearctic') %>% slice(chull(NMDS1, NMDS2))
realm_nmds_afrotropic_polygon = cities_to_realms_nmds %>% filter(core_realm == 'Afrotropic') %>% slice(chull(NMDS1, NMDS2))
realm_nmds_indomalayan_polygon = cities_to_realms_nmds %>% filter(core_realm == 'Indomalayan') %>% slice(chull(NMDS1, NMDS2))
realm_nmds_australasia_polygon = cities_to_realms_nmds %>% filter(core_realm == 'Australasia') %>% slice(chull(NMDS1, NMDS2))
polygon_line_type = 'dashed'
polygon_linewidth = 0.4
with_realms_nmds = function(g) {
g +
geom_polygon(data = realm_nmds_neartic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
geom_polygon(data = realm_nmds_neotropic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
geom_polygon(data = realm_nmds_palearctic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
geom_polygon(data = realm_nmds_afrotropic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
geom_polygon(data = realm_nmds_indomalayan_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
geom_polygon(data = realm_nmds_australasia_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0)
}
analysis_data_latlong_coords = analysis_data[,c('longitude', 'latitude')]
coordinates(analysis_data_latlong_coords) = ~ longitude + latitude
analysis_data_coords_nearneigh <- knearneigh(analysis_data_latlong_coords, longlat = TRUE)
analysis_data_neighbours <- knn2nb(analysis_data_coords_nearneigh)
cities_to_realms_latlong = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv')) %>% left_join(analysis_data) %>% filter(!is.na(latitude))
unique(cities_to_realms_latlong$core_realm)
[1] "Nearctic" "Neotropic" "Palearctic" "Afrotropic" "Indomalayan" "Australasia"
realm_latlong_neartic_polygon = cities_to_realms_latlong %>% filter(core_realm == 'Nearctic') %>% slice(chull(latitude, longitude))
realm_latlong_neotropic_polygon = cities_to_realms_latlong %>% filter(core_realm == 'Neotropic') %>% slice(chull(latitude, longitude))
realm_latlong_palearctic_polygon = cities_to_realms_latlong %>% filter(core_realm == 'Palearctic') %>% slice(chull(latitude, longitude))
realm_latlong_afrotropic_polygon = cities_to_realms_latlong %>% filter(core_realm == 'Afrotropic') %>% slice(chull(latitude, longitude))
realm_latlong_indomalayan_polygon = cities_to_realms_latlong %>% filter(core_realm == 'Indomalayan') %>% slice(chull(latitude, longitude))
realm_latlong_australasia_polygon = cities_to_realms_latlong %>% filter(core_realm == 'Australasia') %>% slice(chull(latitude, longitude))
with_realms_latlong = function(g) {
g +
geom_polygon(data = realm_latlong_neartic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
geom_polygon(data = realm_latlong_neotropic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
geom_polygon(data = realm_latlong_palearctic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
geom_polygon(data = realm_latlong_afrotropic_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
geom_polygon(data = realm_latlong_indomalayan_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0) +
geom_polygon(data = realm_latlong_australasia_polygon, linewidth = polygon_linewidth, linetype = polygon_line_type, alpha = 0)
}
with_realms_latlong(ggplot(analysis_data, aes(x = latitude, y = longitude, colour = mntd_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))
moran.test(analysis_data$mntd_standard, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: analysis_data$mntd_standard
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 7.343, p-value = 0.0000000000001044
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.529792986 -0.003257329 0.005269694
with_realms_latlong(ggplot(analysis_data, aes(x = latitude, y = longitude, colour = mntd_normal)) + geom_point() + normalised_colours_scale + labs(colour = "Normalised response"))
moran.test(analysis_data$mntd_normal, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: analysis_data$mntd_normal
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 7.6496, p-value = 0.00000000000001008
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.552360477 -0.003257329 0.005275586
with_realms_nmds(ggplot(analysis_data, aes(x = NMDS1, y = NMDS2, colour = mntd_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))
moran.test(analysis_data$mntd_standard, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: analysis_data$mntd_standard
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 6.5193, p-value = 0.00000000003531
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.445538119 -0.003257329 0.004739028
with_realms_nmds(ggplot(analysis_data, aes(x = NMDS1, y = NMDS2, colour = mntd_normal)) + geom_point() + standardised_colours_scale + labs(colour = "Normalised response"))
moran.test(analysis_data$mntd_normal, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: analysis_data$mntd_normal
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 6.2392, p-value = 0.0000000002199
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.426489674 -0.003257329 0.004744227
with_realms_latlong(ggplot(analysis_data, aes(x = latitude, y = longitude, colour = beak_width_fdiv_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))
moran.test(analysis_data$beak_width_fdiv_standard, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: analysis_data$beak_width_fdiv_standard
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 7.7208, p-value = 0.000000000000005781
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.557818587 -0.003257329 0.005281064
with_realms_latlong(ggplot(analysis_data, aes(x = latitude, y = longitude, colour = beak_width_fdiv_normal)) + geom_point() + standardised_colours_scale + labs(colour = "Normalised response"))
moran.test(analysis_data$beak_width_fdiv_normal, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: analysis_data$beak_width_fdiv_normal
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 8.3245, p-value < 0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.601742587 -0.003257329 0.005281967
with_realms_nmds(ggplot(analysis_data, aes(x = NMDS1, y = NMDS2, colour = beak_width_fdiv_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))
moran.test(analysis_data$beak_width_fdiv_standard, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: analysis_data$beak_width_fdiv_standard
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 8.0587, p-value = 0.0000000000000003855
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.552096148 -0.003257329 0.004749061
with_realms_nmds(ggplot(analysis_data, aes(x = NMDS1, y = NMDS2, colour = beak_width_fdiv_normal)) + geom_point() + standardised_colours_scale + labs(colour = "Normalised response"))
moran.test(analysis_data$beak_width_fdiv_normal, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: analysis_data$beak_width_fdiv_normal
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 7.5255, p-value = 0.00000000000002625
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.515397518 -0.003257329 0.004749858
with_realms_latlong(ggplot(analysis_data, aes(x = latitude, y = longitude, colour = hwi_fdiv_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))
moran.test(analysis_data$hwi_fdiv_standard, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: analysis_data$hwi_fdiv_standard
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 9.756, p-value < 0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.705943030 -0.003257329 0.005284408
with_realms_latlong(ggplot(analysis_data, aes(x = latitude, y = longitude, colour = hwi_fdiv_normal)) + geom_point() + standardised_colours_scale + labs(colour = "Normalised response"))
moran.test(analysis_data$hwi_fdiv_normal, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: analysis_data$hwi_fdiv_normal
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 10.186, p-value < 0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.737118550 -0.003257329 0.005283050
with_realms_nmds(ggplot(analysis_data, aes(x = NMDS1, y = NMDS2, colour = hwi_fdiv_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))
moran.test(analysis_data$hwi_fdiv_standard, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: analysis_data$hwi_fdiv_standard
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 5.7294, p-value = 0.000000005039
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.391697713 -0.003257329 0.004752012
with_realms_nmds(ggplot(analysis_data, aes(x = NMDS1, y = NMDS2, colour = hwi_fdiv_normal)) + geom_point() + standardised_colours_scale + labs(colour = "Normalised response"))
moran.test(analysis_data$hwi_fdiv_normal, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: analysis_data$hwi_fdiv_normal
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 6.2803, p-value = 0.0000000001689
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.429622019 -0.003257329 0.004750814
with_realms_latlong(ggplot(analysis_data, aes(x = latitude, y = longitude, colour = mass_fdiv_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))
moran.test(analysis_data$mass_fdiv_standard, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: analysis_data$mass_fdiv_standard
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 9.1521, p-value < 0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.661355126 -0.003257329 0.005273416
with_realms_latlong(ggplot(analysis_data, aes(x = latitude, y = longitude, colour = mass_fdiv_normal)) + geom_point() + standardised_colours_scale + labs(colour = "Normalised response"))
moran.test(analysis_data$mass_fdiv_normal, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: analysis_data$mass_fdiv_normal
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 9.008, p-value < 0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.650776207 -0.003257329 0.005271668
with_realms_nmds(ggplot(analysis_data, aes(x = NMDS1, y = NMDS2, colour = mass_fdiv_standard)) + geom_point() + standardised_colours_scale + labs(colour = "Standardised response"))
moran.test(analysis_data$mass_fdiv_standard, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: analysis_data$mass_fdiv_standard
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 6.2114, p-value = 0.0000000002625
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.424488243 -0.003257329 0.004742312
with_realms_nmds(ggplot(analysis_data, aes(x = NMDS1, y = NMDS2, colour = mass_fdiv_normal)) + geom_point() + standardised_colours_scale + labs(colour = "Normalised response"))
moran.test(analysis_data$mass_fdiv_normal, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: analysis_data$mass_fdiv_normal
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 4.6538, p-value = 0.000001629
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.317170749 -0.003257329 0.004740769
all_explanatories = c(
'city_avg_ndvi', 'city_avg_elevation', 'city_avg_temp',
'region_50km_avg_soil_moisture',
'core_realmAfrotropic', 'core_realmAustralasia', 'core_realmIndomalayan', 'core_realmNearctic', 'core_realmNeotropic', 'core_realmPalearctic',
'has_introduced_speciesNo introduced species', 'has_introduced_speciesIntroduced species'
)
all_explanatory_names = factor(
c(
'Avg. NDVI', 'Avg. Elevation', 'Avg. Temp.',
'Avg. Soil Moisture',
'Afrotropic', 'Australasia', 'Indomalayan', 'Nearctic', 'Neotropic', 'Palearctic',
'Introduced Absent', 'Introduced Present'
), ordered = T
)
explanatory_dictionary = data.frame(explanatory = all_explanatories, explanatory_name = all_explanatory_names)
with_explanatory_type_labels = function(p) {
p = p[p$explanatory != '(Intercept)',]
explanatory_levels = all_explanatories[all_explanatories %in% p$explanatory]
p$explanatory <- factor(p$explanatory, levels = explanatory_levels)
p$type <- 'Realm'
p$type[p$explanatory %in% c('city_avg_ndvi', 'city_avg_elevation', 'city_avg_temp')] <- 'City geography'
p$type[p$explanatory %in% c('region_50km_avg_soil_moisture')] <- 'Regional (50 km) geography'
p$type[p$explanatory %in% c('has_introduced_speciesNo introduced species', 'has_introduced_speciesIntroduced species')] <- 'Introduced species'
p
}
with_explanatory_names = function(p) {
p %>% left_join(explanatory_dictionary) %>% arrange(desc(explanatory_name))
}
type_labels = function(p) {
explanatory_levels = all_explanatories[all_explanatories %in% p$explanatory]
p$explanatory <- factor(p$explanatory, levels = explanatory_levels)
p$type <- 'Realm'
p$type[p$explanatory %in% c('city_avg_ndvi', 'city_avg_elevation', 'city_avg_temp', 'city_avg_min_monthly_temp', 'city_avg_max_monthly_temp',
'city_avg_monthly_temp', 'city_avg_rainfall', 'city_avg_max_monthly_rainfall', 'city_avg_min_monthly_rainfall',
'city_avg_soil_moisture', 'city_max_elev', 'city_min_elev', 'city_elev_range')] <- 'City geography'
p$type[p$explanatory %in% c('region_50km_avg_ndvi', 'region_50km_avg_elevation', 'region_50km_avg_soil_moisture', 'region_50km_max_elev',
'region_50km_min_elev', 'region_50km_elev_range')] <- 'Regional (50 km) geography'
p$type[p$explanatory %in% c('region_20km_avg_ndvi', 'region_20km_avg_elevation', 'region_20km_avg_soil_moisture', 'region_20km_max_elev',
'region_20km_min_elev', 'region_20km_elev_range')] <- 'Regional (20 km) geography'
p$type[p$explanatory %in% c('has_introduced_speciesNo introduced species', 'has_introduced_speciesIntroduced species')] <- 'Introduced species'
p
}
explanatory_labels = c(
'has_introduced_species'='Introduced species',
'has_introduced_speciesNo introduced species'='Introduced absent',
'has_introduced_speciesIntroduced species'='Introduced present',
'city_avg_ndvi'='Average NDVI',
'city_avg_elevation'='Average elevation',
'city_avg_temp'='Average temperature',
'city_avg_min_monthly_temp'='Average minimum monthly temperature',
'city_avg_max_monthly_temp'='Average maximum monthly temperature',
'city_avg_monthly_temp'='Average monthly temperature',
'city_avg_rainfall'='Average rainfall',
'city_avg_max_monthly_rainfall'='Average maximum monthly rainfall',
'city_avg_min_monthly_rainfall'='Average minimum monthly rainfall',
'city_avg_soil_moisture'='Average soil moisture',
'city_max_elev'='Maximum elevation',
'city_min_elev'='Minimum elevation',
'city_elev_range'='Elevation range',
'region_20km_avg_ndvi'='Average NDVI',
'region_20km_avg_elevation'='Average elevation',
'region_20km_avg_soil_moisture'='Average soil moisture',
'region_20km_max_elev'='Maximum elevation',
'region_20km_min_elev'='Minimum elevation',
'region_20km_elev_range'='Elevation range',
'region_50km_avg_ndvi'='Average NDVI',
'region_50km_avg_elevation'='Average elevation',
'region_50km_avg_soil_moisture'='Average soil moisture',
'region_50km_max_elev'='Maximum elevation',
'region_50km_min_elev'='Minimum elevation',
'region_50km_elev_range'='Elevation range',
'abs_latitude' = 'Absolute latitude',
'latitude' = 'Latitude',
'longitude' = 'Longitude',
'core_realmAfrotropic' = 'Afrotropical',
'core_realmAustralasia' = 'Austaliasian',
'core_realmIndomalayan' = 'Indomalayan',
'core_realmNearctic' = 'Nearctic',
'core_realmNeotropic' = 'Neotropical',
'core_realmPalearctic' = 'Palearctic',
'core_realmOceania' = 'Oceanical')
create_formula = function(response_var) {
as.formula(paste(response_var, '~ core_realm + city_avg_ndvi + city_avg_elevation + city_avg_temp + region_50km_avg_soil_moisture + has_introduced_species'))
}
geom_map = function(map_sf, title, scale = standardised_colours_scale, colour_label = 'Standardised\nResponse') {
norm_mntd_analysis_geo = ggplot() +
geom_sf(data = world_map, aes(geometry = geometry)) +
map_sf +
scale +
labs(colour = colour_label) +
theme_bw() +
theme(legend.position="bottom")
}
geom_map_std = function(map_sf, title) {
geom_map(map_sf, title)
}
geom_map_nrm = function(map_sf, title) {
geom_map(map_sf, title, normalised_colours_scale, 'Normalised\nResponse')
}
# Taken from MuMIN package
# https://rdrr.io/cran/MuMIn/src/R/averaging.R
# https://rdrr.io/cran/MuMIn/src/R/model.avg.R
.coefarr.avg <-
function(cfarr, weight, revised.var, full, alpha) {
weight <- weight / sum(weight)
nCoef <- dim(cfarr)[3L]
if(full) {
nas <- is.na(cfarr[, 1L, ]) & is.na(cfarr[, 2L, ])
cfarr[, 1L, ][nas] <- cfarr[, 2L, ][nas] <- 0
#cfarr[, 1L:2L, ][is.na(cfarr[, 1L:2L, ])] <- 0
if(!all(is.na(cfarr[, 3L, ])))
cfarr[ ,3L, ][is.na(cfarr[ , 3L, ])] <- Inf
}
avgcoef <- array(dim = c(nCoef, 5L),
dimnames = list(dimnames(cfarr)[[3L]], c("Estimate",
"Std. Error", "Adjusted SE", "Lower CI", "Upper CI")))
for(i in seq_len(nCoef))
avgcoef[i, ] <- par.avg(cfarr[, 1L, i], cfarr[, 2L, i], weight,
df = cfarr[, 3L, i], alpha = alpha, revised.var = revised.var)
avgcoef[is.nan(avgcoef)] <- NA
return(avgcoef)
}
.makecoefmat <- function(cf) {
no.ase <- all(is.na(cf[, 3L]))
z <- abs(cf[, 1L] / cf[, if(no.ase) 2L else 3L])
pval <- 2 * pnorm(z, lower.tail = FALSE)
cbind(cf[, if(no.ase) 1L:2L else 1L:3L, drop = FALSE],
`z value` = z, `Pr(>|z|)` = zapsmall(pval))
}
# Generate model selections using lmer, dredge, and model.avg
# `forumla` : a two-sided linear formula object describing both the fixed-effects and random-effects part of the model
# `data` : the data frame containing the variables from the formula
# `aic_delta` : the AIC delta to use for selecting models in model average
model_average <- function(formula, data, aic_delta = 20) {
model <- lm(
formula,
data=data
)
dredge_result <- dredge(model)
summary(model.avg(dredge_result, subset = delta < aic_delta))
}
# Create a summary data frame containing the selected variables from a model
# `model_sum` : The model summary output from `model_average`
model_summary <- function(model_sum) {
.column_name <- function(postfix) {
postfix
}
# just return the estimate and p value
weight <- model_sum$msTable[, 5L]
coefmat.full <- as.data.frame(.makecoefmat(.coefarr.avg(model_sum$coefArray, weight,
attr(model_sum, "revised.var"), TRUE, 0.05)))
coefmat.subset <-
as.data.frame(.makecoefmat(.coefarr.avg(model_sum$coefArray, weight,
attr(model_sum, "revised.var"), FALSE, 0.05)))
coefmat.subset <- coefmat.subset[-c(1), c(1, 2, 5)]
names(coefmat.subset) <- c(.column_name("estimate"), .column_name("error"), .column_name("p"))
coefmat.subset <- tibble::rownames_to_column(coefmat.subset, "explanatory")
coefmat.subset$model = 'subset'
coefmat.full <- coefmat.full[-c(1), c(1, 2, 5)]
names(coefmat.full) <- c(.column_name("estimate"), .column_name("error"), .column_name("p"))
coefmat.full <- tibble::rownames_to_column(coefmat.full, "explanatory")
coefmat.full$model = 'full'
rbind(coefmat.full, coefmat.subset)
}
plot_dredge_result = function(result_table, mu = 0) {
p = result_table[result_table$model == 'full',]
p = type_labels(p)
ggplot(p, aes(y = explanatory, x = estimate, colour = type)) +
geom_line() +
geom_point() +
geom_errorbar(aes(xmin=estimate-error, xmax=estimate+error), width=.2,
position=position_dodge(0.05)) +
scale_y_discrete(
limits = rev(levels(p$explanatory)),
labels = explanatory_labels) +
scale_colour_manual(
values = c(realm_colour, city_geography_colour, regional_50km_geography_colour, regional_20km_geography_colour, introduced_species_colour),
breaks = c('Realm', 'City geography', 'Regional (50 km) geography', 'Regional (20 km) geography', 'Introduced species')) +
theme_bw() +
geom_vline(xintercept=mu, linetype="dotted") +
guides(colour=guide_legend(title="Predictor type")) + xlab('Difference in response from 0\nhabitat filtering (< 0) and competitive interactions (> 0)\n± Standard Error') + ylab('Predictor') +
theme(legend.justification = "top")
}
gls_method = "ML"
spatial_model = function(formula, correlation) {
gls(
formula,
data = analysis_data,
correlation = correlation,
method = gls_method
)
}
plot_spatial_result = function(model_result) {
model_summary = summary(model_result)
result_table = as.data.frame(model_summary$tTable)
result_table$explanatory = rownames(result_table)
result_table = result_table %>% with_explanatory_type_labels() %>% with_explanatory_names()
ggplot2::ggplot(result_table, ggplot2::aes(y=factor(explanatory_name, level = all_explanatory_names, ordered = T), x=Value, colour = type)) +
ggplot2::geom_line() +
ggplot2::geom_point() +
ggplot2::geom_errorbar(ggplot2::aes(xmin=Value-Std.Error, xmax=Value+Std.Error), width=.2,
position=ggplot2::position_dodge(0.05)) +
ggplot2::theme_bw() +
ggplot2::geom_vline(xintercept=0, linetype="dotted") +
ggplot2::theme(legend.justification = "top") +
ylab('Predictor') +
guides(colour=guide_legend(title="Predictor type")) + xlab('Difference in response from 0\nhabitat filtering (< 0) and competitive interactions (> 0)\n± Standard Error') +
scale_colour_manual(
values = c(realm_colour, city_geography_colour, regional_50km_geography_colour, introduced_species_colour),
breaks = c('Realm', 'City geography', 'Regional (50 km) geography', 'Introduced species')) +
scale_y_discrete(limits = rev(all_explanatory_names[all_explanatory_names %in% result_table$explanatory_name]))
}
AIC(spatial_model(create_formula('mntd_standard'), corLin(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 661.8901
AIC(spatial_model(create_formula('mntd_standard'), corLin(form = ~ latitude + longitude)))
[1] 666.4295
AIC(spatial_model(create_formula('mntd_standard'), corExp(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 653.0099
AIC(spatial_model(create_formula('mntd_standard'), corExp(form = ~ latitude + longitude)))
[1] 656.6118
AIC(spatial_model(create_formula('mntd_standard'), corGaus(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 657.7216
AIC(spatial_model(create_formula('mntd_standard'), corGaus(form = ~ latitude + longitude)))
[1] 662.5602
AIC(spatial_model(create_formula('mntd_standard'), corRatio(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 650.3759
AIC(spatial_model(create_formula('mntd_standard'), corRatio(form = ~ latitude + longitude)))
[1] 654.8339
AIC(spatial_model(create_formula('mntd_standard'), corSpher(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 659.8723
AIC(spatial_model(create_formula('mntd_standard'), corSpher(form = ~ latitude + longitude)))
[1] 664.1159
MNTD: corRatio with NMDS + lat/long
AIC(spatial_model(create_formula('mntd_normal'), corLin(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -81.16061
AIC(spatial_model(create_formula('mntd_normal'), corLin(form = ~ latitude + longitude)))
[1] -77.03513
AIC(spatial_model(create_formula('mntd_normal'), corExp(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -89.28852
AIC(spatial_model(create_formula('mntd_normal'), corExp(form = ~ latitude + longitude)))
[1] -86.17138
AIC(spatial_model(create_formula('mntd_normal'), corGaus(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -86.02288
AIC(spatial_model(create_formula('mntd_normal'), corGaus(form = ~ latitude + longitude)))
[1] -79.92235
AIC(spatial_model(create_formula('mntd_normal'), corRatio(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -92.46759
AIC(spatial_model(create_formula('mntd_normal'), corRatio(form = ~ latitude + longitude)))
[1] -86.60443
AIC(spatial_model(create_formula('mntd_normal'), corSpher(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -86.78558
AIC(spatial_model(create_formula('mntd_normal'), corSpher(form = ~ latitude + longitude)))
[1] -81.47393
AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corLin(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 793.5251
#AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corLin(form = ~ latitude + longitude)))
AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corExp(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 789.282
AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corExp(form = ~ latitude + longitude)))
[1] 790.0085
AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corGaus(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 793.7668
AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corGaus(form = ~ latitude + longitude)))
[1] 793.7656
AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corRatio(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 791.2572
AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corRatio(form = ~ latitude + longitude)))
[1] 791.3899
AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corSpher(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 793.8106
AIC(spatial_model(create_formula('beak_width_fdiv_standard'), corSpher(form = ~ latitude + longitude)))
[1] 793.8106
Beak width: corExp with NMDS + lat/long
AIC(spatial_model(create_formula('beak_width_fdiv_normal'), corLin(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -45.30823
AIC(spatial_model(create_formula('beak_width_fdiv_normal'), corLin(form = ~ latitude + longitude)))
[1] -45.30823
AIC(spatial_model(create_formula('beak_width_fdiv_normal'), corExp(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -51.4648
AIC(spatial_model(create_formula('beak_width_fdiv_normal'), corExp(form = ~ latitude + longitude)))
[1] -50.46273
AIC(spatial_model(create_formula('beak_width_fdiv_normal'), corGaus(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -45.13365
AIC(spatial_model(create_formula('beak_width_fdiv_normal'), corGaus(form = ~ latitude + longitude)))
[1] -45.13552
AIC(spatial_model(create_formula('beak_width_fdiv_normal'), corRatio(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -47.44015
AIC(spatial_model(create_formula('beak_width_fdiv_normal'), corRatio(form = ~ latitude + longitude)))
[1] -47.368
AIC(spatial_model(create_formula('beak_width_fdiv_normal'), corSpher(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -44.97422
AIC(spatial_model(create_formula('beak_width_fdiv_normal'), corSpher(form = ~ latitude + longitude)))
[1] -44.97475
AIC(spatial_model(create_formula('hwi_fdiv_standard'), corLin(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 694.55
AIC(spatial_model(create_formula('hwi_fdiv_standard'), corLin(form = ~ latitude + longitude)))
[1] 699.3174
AIC(spatial_model(create_formula('hwi_fdiv_standard'), corExp(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 666.6631
AIC(spatial_model(create_formula('hwi_fdiv_standard'), corExp(form = ~ latitude + longitude)))
[1] 671.4924
AIC(spatial_model(create_formula('hwi_fdiv_standard'), corGaus(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 704.4854
AIC(spatial_model(create_formula('hwi_fdiv_standard'), corGaus(form = ~ latitude + longitude)))
[1] 709.2567
AIC(spatial_model(create_formula('hwi_fdiv_standard'), corRatio(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 685.5013
AIC(spatial_model(create_formula('hwi_fdiv_standard'), corRatio(form = ~ latitude + longitude)))
[1] 690.9801
AIC(spatial_model(create_formula('hwi_fdiv_standard'), corSpher(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 689.5034
AIC(spatial_model(create_formula('hwi_fdiv_standard'), corSpher(form = ~ latitude + longitude)))
[1] 707.7108
HWI: corExp with NMDS + lat/long
#AIC(spatial_model(create_formula('hwi_fdiv_normal'), corLin(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
#AIC(spatial_model(create_formula('hwi_fdiv_normal'), corLin(form = ~ latitude + longitude)))
AIC(spatial_model(create_formula('hwi_fdiv_normal'), corExp(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -96.05945
AIC(spatial_model(create_formula('hwi_fdiv_normal'), corExp(form = ~ latitude + longitude)))
[1] -91.16342
AIC(spatial_model(create_formula('hwi_fdiv_normal'), corGaus(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -54.57955
AIC(spatial_model(create_formula('hwi_fdiv_normal'), corGaus(form = ~ latitude + longitude)))
[1] -50.52492
AIC(spatial_model(create_formula('hwi_fdiv_normal'), corRatio(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -74.90875
AIC(spatial_model(create_formula('hwi_fdiv_normal'), corRatio(form = ~ latitude + longitude)))
[1] -69.65899
AIC(spatial_model(create_formula('hwi_fdiv_normal'), corSpher(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -53.60829
AIC(spatial_model(create_formula('hwi_fdiv_normal'), corSpher(form = ~ latitude + longitude)))
[1] -61.05117
AIC(spatial_model(create_formula('mass_fdiv_standard'), corLin(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 703.1192
#AIC(spatial_model(create_formula('mass_fdiv_standard'), corLin(form = ~ latitude + longitude)))
AIC(spatial_model(create_formula('mass_fdiv_standard'), corExp(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 753.3038
AIC(spatial_model(create_formula('mass_fdiv_standard'), corExp(form = ~ latitude + longitude)))
[1] 755.7536
AIC(spatial_model(create_formula('mass_fdiv_standard'), corGaus(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 786.6415
AIC(spatial_model(create_formula('mass_fdiv_standard'), corGaus(form = ~ latitude + longitude)))
[1] 785.0663
AIC(spatial_model(create_formula('mass_fdiv_standard'), corRatio(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 766.9612
AIC(spatial_model(create_formula('mass_fdiv_standard'), corRatio(form = ~ latitude + longitude)))
[1] 766.903
AIC(spatial_model(create_formula('mass_fdiv_standard'), corSpher(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] 780.4268
AIC(spatial_model(create_formula('mass_fdiv_standard'), corSpher(form = ~ latitude + longitude)))
[1] 779.8638
Mass: corExp with NMDS + lat/long
AIC(spatial_model(create_formula('mass_fdiv_normal'), corLin(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -98.2492
#AIC(spatial_model(create_formula('mass_fdiv_normal'), corLin(form = ~ latitude + longitude)))
AIC(spatial_model(create_formula('mass_fdiv_normal'), corExp(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -49.779
AIC(spatial_model(create_formula('mass_fdiv_normal'), corExp(form = ~ latitude + longitude)))
[1] -48.0529
AIC(spatial_model(create_formula('mass_fdiv_normal'), corGaus(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -21.84722
AIC(spatial_model(create_formula('mass_fdiv_normal'), corGaus(form = ~ latitude + longitude)))
[1] -23.29214
AIC(spatial_model(create_formula('mass_fdiv_normal'), corRatio(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -35.17179
AIC(spatial_model(create_formula('mass_fdiv_normal'), corRatio(form = ~ latitude + longitude)))
[1] -35.60617
AIC(spatial_model(create_formula('mass_fdiv_normal'), corSpher(form = ~ NMDS1 + NMDS2 + latitude + longitude)))
[1] -27.00252
AIC(spatial_model(create_formula('mass_fdiv_normal'), corSpher(form = ~ latitude + longitude)))
[1] -27.91974
correlation_formula = as.formula('~ NMDS1 + NMDS2 + latitude + longitude')
correlation_function_fdiv = function() {
corExp(form = correlation_formula)
}
correlation_function_mntd = function() {
corRatio(form = correlation_formula)
}
std_mntd_analysis_geo_plot = geom_map_std(geom_sf(data = analysis_data, aes(color = mntd_standard, geometry = geometry)), 'MNTD')
std_mntd_analysis_geo_plot
std_mntd_analysis_data = model_data(analysis_data[!is.na(analysis_data$mntd_standard),], 'mntd_standard')
std_mntd_analysis_formula = create_formula('mntd_standard')
std_mntd_analysis_result = model_average(std_mntd_analysis_formula, std_mntd_analysis_data)
std_mntd_analysis_result_table = model_summary(std_mntd_analysis_result)
std_mntd_analysis_result_table
std_mntd_analysis_pred_plot = plot_dredge_result(std_mntd_analysis_result_table)
std_mntd_analysis_pred_plot
Do the residuals still contain spatial autocorrelation from a fitted lm?
std_mntd_lm = lm(std_mntd_analysis_formula, std_mntd_analysis_data)
moran.test(std_mntd_lm$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: std_mntd_lm$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 3.8425, p-value = 0.0000609
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.275671194 -0.003257329 0.005269438
moran.test(std_mntd_lm$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: std_mntd_lm$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 2.9566, p-value = 0.001555
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.200271681 -0.003257329 0.004738802
std_mntd_spatial_model = spatial_model(std_mntd_analysis_formula, correlation_function_mntd())
moran.test(std_mntd_spatial_model$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: std_mntd_spatial_model$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 4.0485, p-value = 0.00002578
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.290664449 -0.003257329 0.005270844
moran.test(std_mntd_spatial_model$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: std_mntd_spatial_model$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 3.1544, p-value = 0.0008041
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.213916976 -0.003257329 0.004740043
std_mntd_analysis_pred_spatial_plot = plot_spatial_result(std_mntd_spatial_model)
std_mntd_analysis_pred_spatial_plot
nrm_mntd_analysis_geo_plot = geom_map_nrm(geom_sf(data = analysis_data, aes(color = mntd_normal, geometry = geometry)), 'MNTD')
nrm_mntd_analysis_geo_plot
nrm_mntd_analysis_data = model_data(analysis_data[!is.na(analysis_data$mntd_normal),], 'mntd_normal')
nrm_mntd_analysis_formula = create_formula('mntd_normal')
nrm_mntd_analysis_result = model_average(nrm_mntd_analysis_formula, nrm_mntd_analysis_data)
nrm_mntd_analysis_result_table = model_summary(nrm_mntd_analysis_result)
nrm_mntd_analysis_result_table
nrm_mntd_analysis_pred_plot = plot_dredge_result(nrm_mntd_analysis_result_table)
nrm_mntd_analysis_pred_plot
Do the residuals still contain spatial autocorrelation from a fitted lm?
nrm_mntd_lm = lm(nrm_mntd_analysis_formula, nrm_mntd_analysis_data)
moran.test(nrm_mntd_lm$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: nrm_mntd_lm$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 4.2864, p-value = 0.000009079
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.307465451 -0.003257329 0.005254828
moran.test(nrm_mntd_lm$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: nrm_mntd_lm$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 2.0944, p-value = 0.01811
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.140725445 -0.003257329 0.004725908
nrm_mntd_spatial_model = spatial_model(nrm_mntd_analysis_formula, correlation_function_mntd())
moran.test(nrm_mntd_spatial_model$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: nrm_mntd_spatial_model$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 4.7455, p-value = 0.00000104
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.340740843 -0.003257329 0.005254693
moran.test(nrm_mntd_spatial_model$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: nrm_mntd_spatial_model$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 2.3578, p-value = 0.009192
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.158827607 -0.003257329 0.004725789
nrm_mntd_analysis_pred_spatial_plot = plot_spatial_result(nrm_mntd_spatial_model)
nrm_mntd_analysis_pred_spatial_plot
std_gape_fdiv_analysis_geo_plot = geom_map_std(geom_sf(data = analysis_data, aes(color = beak_width_fdiv_standard, geometry = geometry)), 'Beak Width FDiv')
std_gape_fdiv_analysis_geo_plot
std_gape_fdiv_analysis_data = model_data(analysis_data[!is.na(analysis_data$beak_width_fdiv_standard),], 'beak_width_fdiv_standard')
std_gape_fdiv_analysis_formula = create_formula('beak_width_fdiv_standard')
std_gape_fdiv_analysis_result = model_average(std_gape_fdiv_analysis_formula, std_gape_fdiv_analysis_data)
std_gape_fdiv_analysis_result_table = model_summary(std_gape_fdiv_analysis_result)
std_gape_fdiv_analysis_result_table
std_gape_fdiv_analysis_pred_plot = plot_dredge_result(std_gape_fdiv_analysis_result_table)
std_gape_fdiv_analysis_pred_plot
Do the residuals still contain spatial autocorrelation from a fitted lm?
std_gape_fdiv_lm = lm(std_gape_fdiv_analysis_formula, std_gape_fdiv_analysis_data)
moran.test(std_gape_fdiv_lm$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: std_gape_fdiv_lm$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 3.5059, p-value = 0.0002275
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.251017001 -0.003257329 0.005260172
moran.test(std_gape_fdiv_lm$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: std_gape_fdiv_lm$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 1.6582, p-value = 0.04864
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.110794578 -0.003257329 0.004730625
std_gape_fdiv_spatial_model = spatial_model(std_gape_fdiv_analysis_formula, correlation_function_fdiv())
moran.test(std_gape_fdiv_spatial_model$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: std_gape_fdiv_spatial_model$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 3.5695, p-value = 0.0001789
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.255641928 -0.003257329 0.005260840
moran.test(std_gape_fdiv_spatial_model$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: std_gape_fdiv_spatial_model$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 1.868, p-value = 0.03088
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.125233427 -0.003257329 0.004731214
std_gape_fdiv_analysis_pred_spatial_plot = plot_spatial_result(std_gape_fdiv_spatial_model)
std_gape_fdiv_analysis_pred_spatial_plot
nrm_gape_fdiv_analysis_geo_plot = geom_map_nrm(geom_sf(data = analysis_data, aes(color = beak_width_fdiv_normal, geometry = geometry)), 'Beak Width FDiv')
nrm_gape_fdiv_analysis_geo_plot
nrm_gape_fdiv_analysis_data = model_data(analysis_data[!is.na(analysis_data$beak_width_fdiv_normal),], 'beak_width_fdiv_normal')
nrm_gape_fdiv_analysis_formula = create_formula('beak_width_fdiv_normal')
nrm_gape_fdiv_analysis_result = model_average(nrm_gape_fdiv_analysis_formula, nrm_gape_fdiv_analysis_data)
nrm_gape_fdiv_analysis_result_table = model_summary(nrm_gape_fdiv_analysis_result)
nrm_gape_fdiv_analysis_result_table
nrm_gape_fdiv_analysis_pred_plot = plot_dredge_result(nrm_gape_fdiv_analysis_result_table)
nrm_gape_fdiv_analysis_pred_plot
Do the residuals still contain spatial autocorrelation from a fitted lm?
nrm_gape_fdiv_lm = lm(nrm_gape_fdiv_analysis_formula, nrm_gape_fdiv_analysis_data)
moran.test(nrm_gape_fdiv_lm$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: nrm_gape_fdiv_lm$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 4.3357, p-value = 0.000007265
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.311050869 -0.003257329 0.005255227
moran.test(nrm_gape_fdiv_lm$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: nrm_gape_fdiv_lm$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 1.4217, p-value = 0.07756
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.094482089 -0.003257329 0.004726261
nrm_gape_fdiv_spatial_model = spatial_model(nrm_gape_fdiv_analysis_formula, correlation_function_fdiv())
moran.test(nrm_gape_fdiv_spatial_model$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: nrm_gape_fdiv_spatial_model$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 4.4995, p-value = 0.000003406
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.322936781 -0.003257329 0.005255637
moran.test(nrm_gape_fdiv_spatial_model$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: nrm_gape_fdiv_spatial_model$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 1.7347, p-value = 0.0414
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.116002600 -0.003257329 0.004726623
nrm_gape_fdiv_analysis_pred_spatial_plot = plot_spatial_result(nrm_gape_fdiv_spatial_model)
nrm_gape_fdiv_analysis_pred_spatial_plot
std_hwi_fdiv_analysis_geo_plot = geom_map_std(geom_sf(data = analysis_data, aes(color = hwi_fdiv_standard, geometry = geometry)), 'HWI FDiv')
std_hwi_fdiv_analysis_geo_plot
std_hwi_fdiv_analysis_data = model_data(analysis_data[!is.na(analysis_data$hwi_fdiv_standard),], 'hwi_fdiv_standard')
std_hwi_fdiv_analysis_formula = create_formula('hwi_fdiv_standard')
std_hwi_fdiv_analysis_result = model_average(std_hwi_fdiv_analysis_formula, std_hwi_fdiv_analysis_data)
std_hwi_fdiv_analysis_result_table = model_summary(std_hwi_fdiv_analysis_result)
std_hwi_fdiv_analysis_result_table
std_hwi_fdiv_analysis_pred_plot = plot_dredge_result(std_hwi_fdiv_analysis_result_table)
std_hwi_fdiv_analysis_pred_plot
std_hwi_fdiv_lm = lm(std_hwi_fdiv_analysis_formula, std_hwi_fdiv_analysis_data)
moran.test(std_hwi_fdiv_lm$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: std_hwi_fdiv_lm$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 6.1738, p-value = 0.0000000003333
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.444286636 -0.003257329 0.005254904
moran.test(std_hwi_fdiv_lm$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: std_hwi_fdiv_lm$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 3.8009, p-value = 0.00007207
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.258041658 -0.003257329 0.004725975
std_hwi_fdiv_spatial_model = spatial_model(std_hwi_fdiv_analysis_formula, correlation_function_fdiv())
moran.test(std_hwi_fdiv_spatial_model$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: std_hwi_fdiv_spatial_model$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 7.2179, p-value = 0.0000000000002639
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.520272713 -0.003257329 0.005260863
moran.test(std_hwi_fdiv_spatial_model$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: std_hwi_fdiv_spatial_model$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 3.3854, p-value = 0.0003554
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.229604043 -0.003257329 0.004731234
std_hwi_fdiv_analysis_pred_spatial_plot = plot_spatial_result(std_hwi_fdiv_spatial_model)
std_hwi_fdiv_analysis_pred_spatial_plot
nrm_hwi_fdiv_analysis_geo_plot = geom_map_nrm(geom_sf(data = analysis_data, aes(color = hwi_fdiv_normal, geometry = geometry)), 'HWI FDiv')
nrm_hwi_fdiv_analysis_geo_plot
nrm_hwi_fdiv_analysis_data = model_data(analysis_data[!is.na(analysis_data$hwi_fdiv_normal),], 'hwi_fdiv_normal')
nrm_hwi_fdiv_analysis_formula = create_formula('hwi_fdiv_normal')
nrm_hwi_fdiv_analysis_result = model_average(nrm_hwi_fdiv_analysis_formula, nrm_hwi_fdiv_analysis_data)
nrm_hwi_fdiv_analysis_result_table = model_summary(nrm_hwi_fdiv_analysis_result)
nrm_hwi_fdiv_analysis_result_table
nrm_hwi_fdiv_analysis_pred_plot = plot_dredge_result(nrm_hwi_fdiv_analysis_result_table)
nrm_hwi_fdiv_analysis_pred_plot
nrm_hwi_fdiv_lm = lm(nrm_hwi_fdiv_analysis_formula, nrm_hwi_fdiv_analysis_data)
moran.test(nrm_hwi_fdiv_lm$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: nrm_hwi_fdiv_lm$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 7.0294, p-value = 0.000000000001037
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.505183959 -0.003257329 0.005231738
moran.test(nrm_hwi_fdiv_lm$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: nrm_hwi_fdiv_lm$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 3.6271, p-value = 0.0001433
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.245547978 -0.003257329 0.004705532
nrm_hwi_fdiv_spatial_model = spatial_model(nrm_hwi_fdiv_analysis_formula, correlation_function_fdiv())
moran.test(nrm_hwi_fdiv_spatial_model$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: nrm_hwi_fdiv_spatial_model$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 8.1492, p-value < 0.00000000000000022
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.587846115 -0.003257329 0.005261332
moran.test(nrm_hwi_fdiv_spatial_model$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: nrm_hwi_fdiv_spatial_model$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 3.5151, p-value = 0.0002198
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.238537056 -0.003257329 0.004731648
nrm_hwi_fdiv_analysis_pred_spatial_plot = plot_spatial_result(nrm_hwi_fdiv_spatial_model)
nrm_hwi_fdiv_analysis_pred_spatial_plot
std_mass_fdiv_analysis_geo_plot = geom_map_std(geom_sf(data = analysis_data, aes(color = mass_fdiv_standard, geometry = geometry)), 'Mass FDiv')
std_mass_fdiv_analysis_geo_plot
std_mass_fdiv_analysis_data = model_data(analysis_data[!is.na(analysis_data$mass_fdiv_standard),], 'mass_fdiv_standard')
std_mass_fdiv_analysis_formula = create_formula('mass_fdiv_standard')
std_mass_fdiv_analysis_result <- model_average(std_mass_fdiv_analysis_formula, std_mass_fdiv_analysis_data)
std_mass_fdiv_analysis_result_table = model_summary(std_mass_fdiv_analysis_result)
std_mass_fdiv_analysis_result_table
std_mass_fdiv_analysis_pred_plot = plot_dredge_result(std_mass_fdiv_analysis_result_table)
std_mass_fdiv_analysis_pred_plot
std_mass_fdiv_lm = lm(std_mass_fdiv_analysis_formula, std_mass_fdiv_analysis_data)
moran.test(std_mass_fdiv_lm$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: std_mass_fdiv_lm$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 7.4467, p-value = 0.00000000000004785
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.536475629 -0.003257329 0.005253273
moran.test(std_mass_fdiv_lm$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: std_mass_fdiv_lm$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 4.7957, p-value = 0.0000008107
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.326373879 -0.003257329 0.004724537
std_mass_fdiv_spatial_model = spatial_model(std_mass_fdiv_analysis_formula, correlation_function_fdiv())
moran.test(std_mass_fdiv_spatial_model$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: std_mass_fdiv_spatial_model$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 7.9698, p-value = 0.0000000000000007945
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.574414910 -0.003257329 0.005253707
moran.test(std_mass_fdiv_spatial_model$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: std_mass_fdiv_spatial_model$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 5.3014, p-value = 0.00000005745
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.361152294 -0.003257329 0.004724919
std_mass_fdiv_analysis_pred_spatial_plot = plot_spatial_result(std_mass_fdiv_spatial_model)
std_mass_fdiv_analysis_pred_spatial_plot
nrm_mass_fdiv_analysis_geo_plot = geom_map_nrm(geom_sf(data = analysis_data, aes(color = mass_fdiv_normal, geometry = geometry)), 'Mass FDiv')
nrm_mass_fdiv_analysis_geo_plot
nrm_mass_fdiv_analysis_data = model_data(analysis_data[!is.na(analysis_data$mass_fdiv_normal),], 'mass_fdiv_normal')
nrm_mass_fdiv_analysis_formula = create_formula('mass_fdiv_normal')
nrm_mass_fdiv_analysis_result = model_average(nrm_mass_fdiv_analysis_formula, nrm_mass_fdiv_analysis_data)
nrm_mass_fdiv_analysis_result_table = model_summary(nrm_mass_fdiv_analysis_result)
nrm_mass_fdiv_analysis_result_table
nrm_mass_fdiv_analysis_pred_plot = plot_dredge_result(nrm_mass_fdiv_analysis_result_table)
nrm_mass_fdiv_analysis_pred_plot
nrm_mass_fdiv_lm = lm(nrm_mass_fdiv_analysis_formula, nrm_mass_fdiv_analysis_data)
moran.test(nrm_mass_fdiv_lm$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: nrm_mass_fdiv_lm$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 6.9903, p-value = 0.000000000001372
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.503538387 -0.003257329 0.005256267
moran.test(nrm_mass_fdiv_lm$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: nrm_mass_fdiv_lm$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 2.9195, p-value = 0.001753
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.197472475 -0.003257329 0.004727179
nrm_mass_fdiv_spatial_model = spatial_model(nrm_mass_fdiv_analysis_formula, correlation_function_fdiv())
moran.test(nrm_mass_fdiv_spatial_model$residuals, nb2listw(analysis_data_neighbours))
Moran I test under randomisation
data: nrm_mass_fdiv_spatial_model$residuals
weights: nb2listw(analysis_data_neighbours)
Moran I statistic standard deviate = 7.5406, p-value = 0.0000000000000234
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.543587169 -0.003257329 0.005259210
moran.test(nrm_mass_fdiv_spatial_model$residuals, nb2listw(analysis_data_nmds_neighbours))
Moran I test under randomisation
data: nrm_mass_fdiv_spatial_model$residuals
weights: nb2listw(analysis_data_nmds_neighbours)
Moran I statistic standard deviate = 3.2302, p-value = 0.0006185
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.218896280 -0.003257329 0.004729776
nrm_mass_fdiv_analysis_pred_spatial_plot = plot_spatial_result(nrm_mass_fdiv_spatial_model)
nrm_mass_fdiv_analysis_pred_spatial_plot
pred_legend <- ggpubr::get_legend(
# create some space to the left of the legend
std_hwi_fdiv_analysis_pred_plot + theme(legend.box.margin = margin(0, 0, 0, 0)) + guides(colour=guide_legend(ncol=2)) + labs(color = "Predictor type")
)
geo_legend <- ggpubr::get_legend(
# create some space to the left of the legend
std_mass_fdiv_analysis_geo_plot + theme(legend.box.margin = margin(-80, 0, 0, 12), legend.title.position = "top", legend.key.width = unit(10, 'mm')) + labs(color = "Standardised response")
)
legend = plot_grid(
geo_legend,
pred_legend,
nrow = 1
)
legend
plot_grid(
plot_grid(
std_mntd_analysis_geo_plot + theme(legend.position="none"),
std_mntd_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
nrow = 1
) + draw_label("MNTD", size = 16, angle = 90, x = 0.01, y = 0.5),
plot_grid(
std_gape_fdiv_analysis_geo_plot + theme(legend.position="none"),
std_gape_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
nrow = 1
) + draw_label("Beak Width", size = 16, angle = 90, x = 0.01, y = 0.5),
plot_grid(
std_hwi_fdiv_analysis_geo_plot + theme(legend.position="none"),
std_hwi_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
nrow = 1
) + draw_label("HWI", size = 16, angle = 90, x = 0.01, y = 0.5),
plot_grid(
std_mass_fdiv_analysis_geo_plot + theme(legend.position="none"),
std_mass_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
nrow = 1
) + draw_label("Mass", size = 16, angle = 90, x = 0.01, y = 0.5),
legend,
nrow = 5
)
ggsave(filename(FIGURES_OUTPUT_DIR, 'process_response.jpg'), width = 3000, height = 3200, units = 'px')
pred_fig_legend <- ggpubr::get_legend(
# create some space to the left of the legend
std_hwi_fdiv_analysis_pred_plot + theme(legend.box.margin = margin(0, 0, 0, -20)) + guides(colour=guide_legend(ncol=2)) + labs(color = "Predictor type")
)
geo_fig_legend <- ggpubr::get_legend(
# create some space to the left of the legend
std_mass_fdiv_analysis_geo_plot + theme(legend.box.margin = margin(0, 0, 0, 0), legend.title.position = "top", legend.key.width = unit(10, 'mm')) + labs(color = "Standardised response")
)
remove_x_scale = scale_x_continuous(name = '', limits = c(-3, 3))
theme_no_legend = theme(legend.position="none", panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))
fig2 = grid.arrange(
# row 1 - titles
arrangeGrob(grid::textGrob('A) Standardised response by city', x = 0.1, hjust = 0, gp=gpar(fontface="bold"))),
arrangeGrob(grid::textGrob('B) Standardised response predictors', x = 0.1, hjust = 0, gp=gpar(fontface="bold"))),
# row 2
arrangeGrob(
std_mntd_analysis_geo_plot + theme_no_legend,
left = "MNTD"
),
arrangeGrob(
std_mntd_analysis_pred_plot + theme_no_legend + remove_x_scale + ylab('')
),
# row 3
arrangeGrob(std_gape_fdiv_analysis_geo_plot + theme_no_legend, left = "Beak Width"),
arrangeGrob(
std_gape_fdiv_analysis_pred_plot + theme_no_legend + remove_x_scale + ylab('')
),
# row 4
arrangeGrob(std_hwi_fdiv_analysis_geo_plot + theme_no_legend, left = "HWI"),
arrangeGrob(
std_hwi_fdiv_analysis_pred_plot + theme_no_legend + remove_x_scale + ylab('')
),
# row 5
arrangeGrob(std_mass_fdiv_analysis_geo_plot + theme_no_legend, left = "Mass"),
arrangeGrob(
std_mass_fdiv_analysis_pred_plot + theme_no_legend + remove_x_scale + ylab('')
),
# row 6 - legends
arrangeGrob(geo_fig_legend),
arrangeGrob(pred_fig_legend),
heights = c(0.5, 2, 2, 2, 2, 1.25),
nrow = 6
)
jpeg(filename(FIGURES_OUTPUT_DIR, 'figure2.jpg'), width = 183, height = 180, units = 'mm', res = 450)
grid.arrange(
arrangeGrob(fig2),
ncol = 1
)
dev.off()
null device
1
pdf(filename(FIGURES_OUTPUT_DIR, 'figure2.pdf'), width = 12, height = 12, family = "Helvetica")
grid.arrange(
arrangeGrob(fig2),
ncol = 1
)
dev.off()
null device
1
plot_grid(
plot_grid(
ggdraw() +
draw_label(
"Spatial Model",
fontface = 'bold',
x = 0,
hjust = 0
),
ggdraw() +
draw_label(
"Dredge Result",
fontface = 'bold',
x = 0,
hjust = 0
)
),
plot_grid(
std_mntd_analysis_pred_spatial_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
std_mntd_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
nrow = 1
) + draw_label("MNTD", size = 16, angle = 90, x = 0.01, y = 0.5),
plot_grid(
std_gape_fdiv_analysis_pred_spatial_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
std_gape_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
nrow = 1
) + draw_label("Beak Width", size = 16, angle = 90, x = 0.01, y = 0.5),
plot_grid(
std_hwi_fdiv_analysis_pred_spatial_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
std_hwi_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
nrow = 1
) + draw_label("HWI", size = 16, angle = 90, x = 0.01, y = 0.5),
plot_grid(
std_mass_fdiv_analysis_pred_spatial_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
std_mass_fdiv_analysis_pred_plot + theme(legend.position="none") + scale_x_continuous(name = '', limits = c(-3, 3)) + ylab(''),
nrow = 1
) + draw_label("Mass", size = 16, angle = 90, x = 0.01, y = 0.5),
nrow = 5
)
ggsave(filename(FIGURES_OUTPUT_DIR, 'process_response_vs_spatial.jpg'), width = 3000, height = 3200, units = 'px')
ggplot(analysis_data, aes(x = beak_width_fdiv_standard, y = mntd_standard, colour = core_realm)) +
geom_point() +
ylab("MNTD") +
xlab("Beak Width FDiv") +
theme_bw() + labs(color = "Realm")
ggplot(analysis_data, aes(x = hwi_fdiv_standard, y = mntd_standard, colour = core_realm)) +
geom_point() +
ylab("MNTD") +
xlab("HWI FDiv") +
theme_bw() + labs(color = "Realm")
ggplot(analysis_data, aes(x = hwi_fdiv_standard, y = beak_width_fdiv_standard, colour = core_realm)) +
geom_point() +
ylab("Beak Width FDiv") +
xlab("HWI FDiv") +
theme_bw() + labs(color = "Realm")
mntd_fdiv_analysis = analysis_data %>%
dplyr::select(city_id, mntd_standard, hwi_fdiv_standard, beak_width_fdiv_standard, mass_fdiv_standard) %>%
left_join(community_summary) %>%
mutate(urban_pool_perc = urban_pool_size * 100 / regional_pool_size)
mntd_fdiv_analysis
ggpairs(mntd_fdiv_analysis %>% dplyr::select(mntd_standard, hwi_fdiv_standard, beak_width_fdiv_standard, mass_fdiv_standard, regional_pool_size, urban_pool_size, urban_pool_perc), columnLabels = c('MNTD', 'HWI FD', 'Bk FD', 'Mss FD', 'Region Rich.', 'Urban Rich.', '% Urban'))
ggsave(filename(FIGURES_OUTPUT_DIR, 'appendix_standarised_correlation.jpg'))